zDif <- function(beta1, beta2, se1, se2){
	dif <- beta1 - beta2
	se <- sqrt(se1^2 + se2^2)
	ztest <- dif/se
	out <- list(difference = dif, SE = se, zscore = ztest)
	return(out)
}

polyH1 <- function(yearElection, p){
	out.cent.men.opt <- rdrobust(y = subset(d1, codeyear %in% centParties & ANO_ELEICAO %in% yearElection)$avg_male_DF
			, x = subset(d1, codeyear %in% centParties  & ANO_ELEICAO %in% yearElection)$run_var
			, p = p
			, c = 0	
			, cluster = subset(d1, codeyear %in% centParties & ANO_ELEICAO %in% yearElection)$codeyear
			, all = TRUE)
	out.cent.women.opt <- rdrobust(y = subset(d1, codeyear %in% centParties & ANO_ELEICAO %in% yearElection)$avg_fem_DF
			, x = subset(d1, codeyear %in% centParties  & ANO_ELEICAO %in% yearElection)$run_var
			, p = p
			, c = 0	
			, cluster = subset(d1, codeyear %in% centParties & ANO_ELEICAO %in% yearElection)$codeyear
			, all = TRUE)
	zOut <- zDif(beta1 = out.cent.men.opt$coef[3,], beta2 = out.cent.women.opt$coef[3,]
		, se1 = out.cent.men.opt$se[3,], se2 = out.cent.women.opt$se[3,])

	pvalueDif <- ifelse(zOut$zscore > 0, (1-pnorm(zOut$zscore))* 2, pnorm(zOut$zscore) * 2)
	# To return
	return(list(men = summary(out.cent.men.opt), women = summary(out.cent.women.opt), ztest = zOut, pv_zOut = pvalueDif))
}

plotGap <- function(betas, lb, ub, yub, ylb){
	par(mar = c(5, 7, 2, 2))
	plot(x = seq(0.5, 100, 0.5), y = betas, type = 'l', 
		ylim = c(ylb, yub), xlab = 'Bandwidth',
		ylab = 'Gender Gap',
		cex.lab = 1.75, axes = FALSE)
	lines(x = seq(0.5, 100, 0.5), y = lb, lty = 3, lwd = 1.25)
	lines(x = seq(0.5, 100, 0.5), y = ub, lty = 3, lwd = 1.25)
	abline(h = 0, lty = 2, lwd = 0.5)
	axis(1, cex.axis = 1.5)
	axis(2, cex.axis = 1.5)
}


plotTtest <- function(betas, lb, ub, dvName, yub, ylb){
	par(mar = c(5, 7, 2, 2))
	plot(x = seq(0.5, 100, 0.5), y = betas, type = 'l', 
		ylim = c(ylb, yub), xlab = 'Bandwidth',
		ylab = paste0('Avg. Vote Share, ', dvName, ' (t)'),
		cex.lab = 1.75, axes = FALSE)
	lines(x = seq(0.5, 100, 0.5), y = lb, lty = 3, lwd = 1.25)
	lines(x = seq(0.5, 100, 0.5), y = ub, lty = 3, lwd = 1.25)
	abline(h = 0, lty = 2, lwd = 0.5)
	axis(1, cex.axis = 1.5)
	axis(2, cex.axis = 1.5)
}

genSensitivityH1 <- function(yearElection, yearPlot, file){
	pdf(file)
	outMen <- sensitivity(mydata = subset(d1, codeyear %in% centParties & ANO_ELEICAO == yearElection), 
		DV = 'avg_male_DF', 
		ylb = -4, 
		yub = 4, 
		dvName = 'Federal Deputy - Men', 
		clusterVar = subset(d1, codeyear %in% centParties & ANO_ELEICAO == yearElection)$codeyear)
	outWomen <- sensitivity(mydata = subset(d1, codeyear %in% centParties & ANO_ELEICAO == yearElection), 
		DV = 'avg_fem_DF', 
		ylb = -4, 
		yub = 4, 
		dvName = 'Federal Deputy - Women', 
		clusterVar = subset(d1, codeyear %in% centParties & ANO_ELEICAO == yearElection)$codeyear)
	zscores <- zDif(beta1 = outMen$beta, 
		beta2 = outWomen$beta, 
		se1 = outMen$se, 
		se2 = outWomen$se)$zscore
	difference <- zDif(beta1 = outMen$beta, 
		beta2 = outWomen$beta, 
		se1 = outMen$se, 
		se2 = outWomen$se)$difference
	pvalues <- ifelse(zscores > 0, (1 - pnorm(zscores)) * 2, pnorm(zscores) * 2)

	par(mar = c(5, 5, 2, 2))
	plot(x = outWomen$h[difference>0], 
		y = pvalues[difference>0], 
		pch = 20, 
		cex = 1, 
		xlim = c(0, 100),
		ylim = c(0, 1),
		xlab = 'Bandwidth', 
		ylab = 'p-value',
		axes = FALSE, 
		cex.lab = 2)
	points(x = outWomen$h[difference<=0], 
		y = pvalues[difference<=0], col = 'grey40', pch = 18) 
	segments(x0 = -5, x1 = 105, y0 = 0.05, y1 = 0.05, lty = 3)
	segments(x0 = -5, x1 = 105, y0 = 0.10, y1 = 0.10, lty = 2)
	axis(1, at = c(0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100), cex.axis = 1.5)
	axis(2, cex.axis = 1.5)
	legend('topright', 
		legend = c(expression("Men" > "Women"), expression("Men" <= "Women")), 
		col = c('black', 'grey40'),
		pch = c(20, 18), 
		bty = 'n',
		cex = 1.5)
	dev.off()
}


# Function to run sensitivity
sensitivity <- function(mydata, DV, ylb, yub, dvName, clusterVar, run_var = 'run_var'){
	betas <- c()
	lb <- list()
	ub <- list()
	se <- list()
	h <- seq(0.5, 100, by = 0.5)
	for(i in 1:length(h)){
		out <- rdrobust(y = mydata[, DV]
				, x = mydata[, run_var]
				, h = h[i]	
				, cluster = clusterVar		
				, c = 0
				, level = 95
				, all = TRUE)	
		betas[i] <- out$Estimate[,'tau.bc']
		se[i] <- out$Estimate[, 'se.rb']
		lb[i] <- out$ci[3, 1]
		ub[i] <- out$ci[3, 2]
	}
	out <- data.frame(h = h, beta = unlist(betas), se = unlist(se))
	# Create Plot
	par(mar = c(5, 7, 2, 2))
	plot(x = h, y = betas, type = 'l', 
		ylim = c(ylb, yub), xlab = 'Bandwidth',
		ylab = paste0('Avg. Vote Share, ', dvName, ' (t)'),
		cex.lab = 1.75, axes = FALSE)
	lines(x = h, y = unlist(lb), lty = 3, lwd = 1.25)
	lines(x = h, y = unlist(ub), lty = 3, lwd = 1.25)
	abline(h = 0, lty = 2, lwd = 0.5)
	axis(1, cex.axis = 1.5)
	axis(2, cex.axis = 1.5)
	# abline(v = out$bws[1, 1], col = 'red')
	return(out)
}


# Function to generate plots for the paper
genPlot <- function(mydata, DV, dvName, vs = TRUE, fem = FALSE, fed = TRUE, ci = 95, clusterVar, time = ' (t)'){
	d <- mydata

	if(fem){
		mydata <- mydata[mydata$mayorFem == TRUE, ]
	}else{
		mydata <- mydata
	}

	if(fed){
		nonZero <- 'nonZeroDF'
	}else{
		nonZero <- 'nonZeroDE'
	}


	if(vs){
		vs <- 'Avg. Vote Share'
		ylim <- c(0, 5)
		t1 <- 4
		t2 <- 3.65
	}else{
		vs <- 'Total Vote Share'
		ylim <- c(0, 40)
		t1 <- 35
		t2 <- 33		
	}

	m1 <- rdrobust(y = mydata[, DV]
			, x = mydata$run_var
			, p = 1
			, c = 0		
			, level = ci 
			, cluster = clusterVar
			, all = TRUE)
	p1 <- rdplot(y = mydata[, DV]
			, x = mydata$run_var
			)

	par(mar = c(5, 5, 2, 2))
	plot(x = p1$vars_bins$rdplot_mean_x[p1$vars_bins$rdplot_mean_x > -55 & p1$vars_bins$rdplot_mean_x < 55],
		y = p1$vars_bins$rdplot_mean_y[p1$vars_bins$rdplot_mean_x > -55 & p1$vars_bins$rdplot_mean_x < 55],
		ylim = ylim,
		xlim = c(-50, 50), 
		pch = 20, 
		col = 'grey',
		xlab = 'Difference in Vote Share, Municipal Election (t-2)',
		ylab = paste0(vs, ', ', dvName, time),
		axes = FALSE,
		cex.lab = 1.5
		)
	lines(x = p1$vars_poly$rdplot_x[p1$vars_poly$rdplot_x < 0]
		, y = p1$vars_poly$rdplot_y[p1$vars_poly$rdplot_x < 0], lwd = 1.5)
	lines(x = p1$vars_poly$rdplot_x[p1$vars_poly$rdplot_x > 0]
		, y = p1$vars_poly$rdplot_y[p1$vars_poly$rdplot_x > 0], lwd = 1.5)
	axis(1, at = c(seq(-50, 50, by = 25)), cex.axis = 1.25)
	axis(2, cex.axis = 1.25)



	text(x = -20, y = t1, paste0('Estimated Effect = ', round(m1$Estimate[1, 'tau.bc'], 3)), cex = 1.75)
	text(x = -20, y = t2, paste0('p-value = ', ifelse(nchar(round(m1$pv[3, 1], 3)) == 1, paste0(round(m1$pv[3, 1], 3), '.000'), round(m1$pv[3, 1], 3))), cex = 1.75)
}

# Function to run models by party
byParty <- function(party, DV, fed = TRUE){

	if(fed){
		nonZero <- 'nonZeroDF'
	}else{
		nonZero <- 'nonZeroDE'
	}


	d1 <- subset(d, NUMERO_PARTIDO == party)

	m1 <- rdrobust(y = d1[, DV][d1[, nonZero] == 2]
			, x = d1$run_var[d1[, nonZero] == 2]
			, p = 1
			, c = 0		
			, level = 95
			, all = TRUE)
	return(data.frame(party = party
		, N = nrow(d1[d1[, nonZero] == 2,])
		, beta = m1$coef[3, 1]
		, se = m1$se[3, 1]
		, lb = m1$ci[3, 1]
		, ub = m1$ci[3, 2])
	)
}

# Function to run models by year
byYear <- function(year, DV, fed = TRUE){

	if(fed){
		nonZero <- 'nonZeroDF'
	}else{
		nonZero <- 'nonZeroDE'
	}


	d1 <- subset(d, ANO_ELEICAO == year)

	m1 <- rdrobust(y = d1[, DV][d1[, nonZero] == 2]
			, x = d1$run_var[d1[, nonZero] == 2]
			, p = 1
			, c = 0		
			, level = 95
			, all = TRUE)
	return(data.frame(year = year
		, N = nrow(d1[d1[, nonZero] == 2,])
		, beta = m1$coef[3, 1]
		, se = m1$se[3, 1]
		, lb = m1$ci[3, 1]
		, ub = m1$ci[3, 2])
	)
}

genResults <- function(mydata, DV, p = 1, ci = 95, clusterVar){
	d1 <- mydata

	m1 <- rdrobust(y = d1[, DV]
			, x = d1$run_var
			, p = p
			, c = 0		
			, level = ci
			, masspoints = FALSE
			, cluster = clusterVar
			, all = TRUE)

	return(m1)

}


# Start the plot:
byYearH1 <- function(betas_men, betas_women, lb_men, ub_men, lb_women, ub_women, zout, ylims = c(-3, 3), betaPos = 0){

	par(mar = c(5, 5, 2, 2))
	plot(x = seq(2002-1, 2018-1, 4), y = betas_men, axes = FALSE, 
		pch = 20, cex = 2, xlim = c(2000, 2020), ylim = ylims,
		ylab = "RD Effect of Co-Partisan Mayor",
		xlab = "Federal Election", 
		cex.lab = 2)
	sapply(1:5, function(i) lines(x = c(seq(2002-1, 2018-1, 4)[i], seq(2002-1, 2018-1, 4)[i]), 
		y = c(lb_men[i], ub_men[i])))
	sapply(1:5, function(i) lines(x = c(seq(2002+1, 2018+1, 4)[i], seq(2002+1, 2018+1, 4)[i]), 
		y = c(lb_women[i], ub_women[i]), col = 'gray46'))
	axis(1, seq(2002, 2018, 4), at = seq(2002, 2018, 4), tick = FALSE, cex.axis = 2)
	axis(2, tick = T, cex.axis = 2, las = 2)
	segments(x0 = 1997, x1 = 2020.5, y0 = 0, y1 = 0, lty = 3)

	# Add differences
	# # 1998
	# segments(x0 = 1997, x1 = 1998.5, y0 = betas_men[1], y1 = betas_men[1], lty = 2)
	# segments(x0 = 1998.5, x1 = 1998.5, y0 = betas_men[1], y1 = betas_women[1], lty = 2)
	# segments(x0 = 1998.5, x1 = 1999, y0 = betas_women[1], y1 = betas_women[1], lty = 2)
	# 2002
	segments(x0 = 2001, x1 = 2002.5, y0 = betas_men[1], y1 = betas_men[1], lty = 2)
	segments(x0 = 2002.5, x1 = 2002.5, y0 = betas_men[1], y1 = betas_women[1], lty = 2)
	segments(x0 = 2002.5, x1 = 2003, y0 = betas_women[1], y1 = betas_women[1], lty = 2)
	# 2006
	segments(x0 = 2005, x1 = 2006.5, y0 = betas_men[2], y1 = betas_men[2], lty = 2)
	segments(x0 = 2006.5, x1 = 2006.5, y0 = betas_men[2], y1 = betas_women[2], lty = 2)
	segments(x0 = 2006.5, x1 = 2007, y0 = betas_women[2], y1 = betas_women[2], lty = 2)
	# 2010
	segments(x0 = 2009, x1 = 2010.5, y0 = betas_men[3], y1 = betas_men[3], lty = 2)
	segments(x0 = 2010.5, x1 = 2010.5, y0 = betas_men[3], y1 = betas_women[3], lty = 2)
	segments(x0 = 2010.5, x1 = 2011, y0 = betas_women[3], y1 = betas_women[3], lty = 2)
	# 2014
	segments(x0 = 2013, x1 = 2014.5, y0 = betas_men[4], y1 = betas_men[4], lty = 2)
	segments(x0 = 2014.5, x1 = 2014.5, y0 = betas_men[4], y1 = betas_women[4], lty = 2)
	segments(x0 = 2014.5, x1 = 2015, y0 = betas_women[4], y1 = betas_women[4], lty = 2)
	# 2018
	segments(x0 = 2017, x1 = 2018.5, y0 = betas_men[5], y1 = betas_men[5], lty = 2)
	segments(x0 = 2018.5, x1 = 2018.5, y0 = betas_men[5], y1 = betas_women[5], lty = 2)
	segments(x0 = 2018.5, x1 = 2019, y0 = betas_women[5], y1 = betas_women[5], lty = 2)

	# Add points for women's coefficients
	points(x = seq(2002+1, 2018+1, 4), y = betas_women, pch = 17, col = 'gray46', cex = 1.5)

	# Add information about difference
	# text(round(zout[[1]]$difference, 3), y = betas_men[1] + 0.5, x = 1998, cex = 1.5)
	# text(paste0('(', round(ifelse(zout[[1]]$zscore > 0, (1 - pnorm(zout[[1]]$zscore)) * 2, pnorm(zout[[1]]$zscore) * 2), 3), ')')
	# , y = betas_men[1] + 0.2, x = 1998, cex = 1.5)

	text(round(zout[[1]]$difference, 3), y = betas_men[1] + 0.75+ betaPos, x = 2002, cex = 1.5)
	text(paste0('(', round(ifelse(zout[[2]]$zscore > 0, (1 - pnorm(zout[[1]]$zscore)) * 2, pnorm(zout[[1]]$zscore) * 2), 3), ')')
	, y = betas_men[1] + 0.45, x = 2002, cex = 1.5)

	text(round(zout[[2]]$difference, 3), y = betas_men[2] + 0.5+ betaPos, x = 2006, cex = 1.5)
	text(paste0('(', round(ifelse(zout[[3]]$zscore > 0, (1 - pnorm(zout[[2]]$zscore)) * 2, pnorm(zout[[2]]$zscore) * 2), 3), ')')
	, y = betas_men[2] + 0.2, x = 2006, cex = 1.5)

	text(round(zout[[3]]$difference, 3), y = betas_men[3] + 0.5+ betaPos, x = 2010, cex = 1.5)
	text(paste0('(', round(ifelse(zout[[3]]$zscore > 0, (1 - pnorm(zout[[3]]$zscore)) * 2, pnorm(zout[[3]]$zscore) * 2), 3), ')')
	, y = betas_men[3] + 0.2, x = 2010, cex = 1.5)

	text(round(zout[[4]]$difference, 3), y = betas_men[4] + 0.5+ betaPos, x = 2014, cex = 1.5)
	text(paste0('(', round(ifelse(zout[[4]]$zscore > 0, (1 - pnorm(zout[[4]]$zscore)) * 2, pnorm(zout[[4]]$zscore) * 2), 3), ')')
	, y = betas_men[4] + 0.2, x = 2014, cex = 1.5)

	text(round(zout[[5]]$difference, 3), y = betas_men[5] + 0.5+ betaPos, x = 2018, cex = 1.5)
	text(paste0('(', round(ifelse(zout[[5]]$zscore > 0, (1 - pnorm(zout[[5]]$zscore)) * 2, pnorm(zout[[5]]$zscore) * 2), 3), ')')
	, y = betas_men[5] + 0.2, x = 2018, cex = 1.5)

	# Add legend
	legend('bottomleft', 
		legend = c('Women', 'Men'), 
		col = c('gray46', 'black'), 
		pch = c(17, 20), 
		lty = c(1, 1),
		bty = 'n',
		cex = 1.5)
}


linerEstimator <- function(yearElection){
	# Objects to store results
	coefMenBef08 <- coefWomenBef08 <- coefDifBef08 <- c()
	seMenBef08 <- seWomenBef08 <- seDifBef08 <- c()
	i <- 1
	for(bw in seq(0.5, 100, 0.5)){
		dataSmall <- subset(d2, abs(run_var) < bw & ANO_ELEICAO %in% yearElection)
		m1 <- lm(avg_DF ~ winner * men + run_var, data = dataSmall)
		# Before 2008
		coefWomenBef08[i] <- coeftest(m1, vcovCL(m1, cluster = dataSmall$codeyear))['winner', 1]
		coefMenBef08[i] <- coeftest(m1, vcovCL(m1, cluster = dataSmall$codeyear))['winner', 1] + coeftest(m1, vcovCL(m1, cluster = dataSmall$codeyear))['winner:men', 1]
		coefDifBef08[i] <- coeftest(m1, vcovCL(m1, cluster = dataSmall$codeyear))['winner:men', 1]
		seWomenBef08[i] <- sqrt(diag(vcovCL(m1, cluster = dataSmall$codeyear)))['winner']
		seMenBef08[i] <- sqrt(vcovCL(m1, cluster = dataSmall$codeyear)['winner', 'winner'] + vcovCL(m1, cluster = dataSmall$codeyear)['winner:men', 'winner:men'] + 2 * vcovCL(m1, cluster = dataSmall$codeyear)['winner', 'winner:men'])
		seDifBef08[i] <- sqrt(diag(vcovCL(m1, cluster = dataSmall$codeyear)))['winner:men']
		# iterate
		i <- i + 1
	}

	lbMenBefore <- seMenBef08 * qnorm(0.025) + coefMenBef08
	ubMenBefore <- seMenBef08 * qnorm(0.975) + coefMenBef08
	lbWomenBefore <- seWomenBef08 * qnorm(0.025) + coefWomenBef08
	ubWomenBefore <- seWomenBef08 * qnorm(0.975) + coefWomenBef08
	lbDifBefore <- seDifBef08 * qnorm(0.025) + coefDifBef08
	ubDifBefore <- seDifBef08 * qnorm(0.975) + coefDifBef08

	return(list(betaMen = coefMenBef08, lbMen = lbMenBefore, ubMen = ubMenBefore,
		betaWomen = coefWomenBef08, lbWomen = lbWomenBefore, ubWomen = ubWomenBefore,
		betaDif = coefDifBef08, lbDif = lbDifBefore, ubDif = ubDifBefore))
}


genSensitivityH2 <- function(yearElection, run_var, period, file){
	pdf(file)
	outMen <- sensitivity(mydata = subset(d2, ANO_ELEICAO %in% yearElection), 
		DV = 'avg_male_DF', 
		ylb = -10, 
		yub = 10, 
		run_var = run_var,
		dvName = 'Federal Deputy - Men', 
		clusterVar = subset(d2, ANO_ELEICAO %in% yearElection)$codeyear)
	outWomen <- sensitivity(mydata = subset(d2, ANO_ELEICAO %in%  yearElection), 
		DV = 'avg_fem_DF', 
		ylb = -10, 
		yub = 10, 
		run_var = run_var,
		dvName = 'Federal Deputy - Women', 
		clusterVar = subset(d2, ANO_ELEICAO %in% yearElection)$codeyear)
	zscores <- zDif(beta1 = outMen$beta, 
		beta2 = outWomen$beta, 
		se1 = outMen$se, 
		se2 = outWomen$se)$zscore
	difference <- zDif(beta1 = outMen$beta, 
		beta2 = outWomen$beta, 
		se1 = outMen$se, 
		se2 = outWomen$se)$difference
	pvalues <- ifelse(zscores > 0, (1 - pnorm(zscores)) * 2, pnorm(zscores) * 2)

	par(mar = c(5, 5, 2, 2))
	plot(x = outWomen$h[difference>0], 
		y = pvalues[difference>0], 
		pch = 20, 
		cex = 1, 
		xlim = c(0, 100),
		ylim = c(0, 1),
		xlab = 'Bandwidth', 
		ylab = 'p-value',
		axes = FALSE, 
		cex.lab = 2)
	points(x = outWomen$h[difference<=0], 
		y = pvalues[difference<=0], col = 'grey40', pch = 18) 
	segments(x0 = -5, x1 = 105, y0 = 0.05, y1 = 0.05, lty = 3)
	segments(x0 = -5, x1 = 105, y0 = 0.10, y1 = 0.10, lty = 2)
	axis(1, at = c(0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100), cex.axis = 1.5)
	axis(2, cex.axis = 1.5)
	legend('topright', 
		legend = c(expression("Men" > "Women"), expression("Men" <= "Women")), 
		col = c('black', 'grey40'),
		pch = c(20, 18), 
		bty = 'n',
		cex = 1.5)
	dev.off()
}


polyH2 <- function(yearElection, p, run_var){
	out.cent.men.opt <- rdrobust(y = subset(d2, ANO_ELEICAO %in% yearElection)$avg_male_DF
			, x = subset(d2, ANO_ELEICAO %in% yearElection)[, run_var]
			, p = p
			, c = 0	
			, cluster = subset(d2, ANO_ELEICAO %in% yearElection)$codeyear
			, all = TRUE)
	out.cent.women.opt <- rdrobust(y = subset(d2, ANO_ELEICAO %in% yearElection)$avg_fem_DF
			, x = subset(d2, ANO_ELEICAO %in% yearElection)[, run_var]
			, p = p
			, c = 0	
			, cluster = subset(d2, ANO_ELEICAO %in% yearElection)$codeyear
			, all = TRUE)
	zOut <- zDif(beta1 = out.cent.men.opt$coef[3,], beta2 = out.cent.women.opt$coef[3,]
		, se1 = out.cent.men.opt$se[3,], se2 = out.cent.women.opt$se[3,])

	pvalueDif <- ifelse(zOut$zscore > 0, (1-pnorm(zOut$zscore))* 2, pnorm(zOut$zscore) * 2)
	# To return
	return(list(men = summary(out.cent.men.opt), women = summary(out.cent.women.opt), ztest = zOut, pv_zOut = pvalueDif))
}

linearEstimatorH2 <- function(yearElection, DV, run_var){
	m <- as.formula(paste0(DV, "~ winner +", run_var))
	# Objects to store results
	coefBef  <- c()
	seBef  <- c()
	i <- 1
	for(bw in seq(0.5, 100, 0.5)){
		dataSmall <- d2[d2[, run_var] < bw & d2$ANO_ELEICAO %in% yearElection, ]
		m1 <- lm(m, data = dataSmall)
		# Before 2008
		coefBef[i] <- coeftest(m1, vcovCL(m1, cluster = dataSmall$codeyear))['winner', 1]
		seBef[i] <- sqrt(vcovCL(m1, cluster = dataSmall$codeyear)['winner', 'winner']) 
		# iterate
		i <- i + 1
	}

	lbBefore <- seBef * qnorm(0.025) + coefBef
	ubBefore <- seBef * qnorm(0.975) + coefBef

	return(list(beta = coefBef, se = seBef, lb = lbBefore, ub = ubBefore))
}


plotDifpvalue <- function(difference, pvalues){
	par(mar = c(5, 5, 2, 2))
	plot(x = seq(0.5, 100, by = 0.5)[difference>0], 
		y = pvalues[difference>0], 
		pch = 20, 
		cex = 1, 
		xlim = c(0, 100),
		ylim = c(0, 1),
		xlab = 'Bandwidth', 
		ylab = 'p-value',
		axes = FALSE, 
		cex.lab = 2)
	points(x = seq(0.5, 100, by = 0.5)[difference<=0], 
		y = pvalues[difference<=0], col = 'grey40', pch = 18) 
	segments(x0 = -5, x1 = 105, y0 = 0.05, y1 = 0.05, lty = 3)
	segments(x0 = -5, x1 = 105, y0 = 0.10, y1 = 0.10, lty = 2)
	axis(1, at = c(0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100), cex.axis = 1.5)
	axis(2, cex.axis = 1.5)
	legend('topright', 
		legend = c(expression("Men" > "Women"), expression("Men" <= "Women")), 
		col = c('black', 'grey40'),
		pch = c(20, 18), 
		bty = 'n',
		cex = 1.5)
}

polyGap <- function(yearElection, p){
	out.cent.opt <- rdrobust(y = subset(d1, codeyear %in% centParties & ANO_ELEICAO == yearElection)$gap
			, x = subset(d1, codeyear %in% centParties  & ANO_ELEICAO == yearElection)$run_var
			, p = p
			, c = 0	
			, cluster = subset(d1, codeyear %in% centParties & ANO_ELEICAO == yearElection)$codeyear
			, all = TRUE)
	# To return
	return(summary(out.cent.opt))
}


linearEstimatorGap <- function(yearElection){
	# Objects to store results
	coefAll  <- c() 
	seAll  <- c()
	i <- 1
	for(bw in seq(0.5, 100, 0.5)){
		dataSmall <- subset(d1, abs(run_var) < bw & ANO_ELEICAO == yearElection)
		m1 <- lm(gap ~ winner + run_var, data = dataSmall)
		coefAll[i] <- coeftest(m1, vcovCL(m1, cluster = dataSmall$codeyear))['winner', 1]
		seAll[i] <- sqrt(diag(vcovCL(m1, cluster = dataSmall$codeyear)))['winner']
		# iterate
		i <- i + 1
	}

	lbAll <- seAll * qnorm(0.025) + coefAll
	ubAll <- seAll * qnorm(0.975) + coefAll

	return(list(beta= coefAll, lb = lbAll, ub = ubAll))
}



# Start the plot:
byPeriodH1 <- function(betas_men, betas_women, lb_men, ub_men, lb_women, ub_women, zout, ylims = c(-3, 3), betaPos = 0){

	par(mar = c(5, 5, 2, 2))
	plot(x = seq(1-0.1, 3-0.1, 1), y = betas_men, axes = FALSE, 
		pch = 20, cex = 2, xlim = c(0.75, 3.25), ylim = ylims,
		ylab = "RD Effect of Co-Partisan Mayor",
		xlab = "Electoral Rule Regime", 
		cex.lab = 2)	
	sapply(1:3, function(i) lines(x = c(seq(1-0.1, 3-0.1, 1)[i], seq(1-0.1, 3-0.1, 1)[i]), 
		y = c(lb_men[i], ub_men[i])))
	sapply(1:3, function(i) lines(x = c(seq(1+0.1, 3+0.1, 1)[i], seq(1+0.1, 3+0.1, 1)[i]), 
		y = c(lb_women[i], ub_women[i]), col = 'gray46'))
	axis(1, at = seq(1, 3, 1), label = c('Weak', 'Strong', 'Moderate'), tick = FALSE, cex.axis = 2)
	axis(2, tick = T, cex.axis = 2, las = 2)
	segments(x0 = 0, x1 = 3.5, y0 = 0, y1 = 0, lty = 3)


	# Add differences
	# Weak
	segments(x0 = 0.9, x1 = 1.05, y0 = betas_men[1], y1 = betas_men[1], lty = 2)
	segments(x0 = 1.05, x1 = 1.05, y0 = betas_men[1], y1 = betas_women[1], lty = 2)
	segments(x0 = 1.05, x1 = 1.1, y0 = betas_women[1], y1 = betas_women[1], lty = 2)
	# Strong
	segments(x0 = 1.9, x1 = 2.05, y0 = betas_men[2], y1 = betas_men[2], lty = 2)
	segments(x0 = 2.05, x1 = 2.05, y0 = betas_men[2], y1 = betas_women[2], lty = 2)
	segments(x0 = 2.05, x1 = 2.1, y0 = betas_women[2], y1 = betas_women[2], lty = 2)
	# Moderate
	segments(x0 = 2.9, x1 = 3.05, y0 = betas_men[3], y1 = betas_men[3], lty = 2)
	segments(x0 = 3.05, x1 = 3.05, y0 = betas_men[3], y1 = betas_women[3], lty = 2)
	segments(x0 = 3.05, x1 = 3.1, y0 = betas_women[3], y1 = betas_women[3], lty = 2)

	# Add points for women's coefficients
	points(x = seq(1+0.1, 3+0.1, 1), y = betas_women, pch = 17, col = 'gray46', cex = 1.5)

	# Add information about difference
	text(round(zout[[1]]$difference, 3), y = betas_men[1] + 0.7, x = 1, cex = 1.5)
	text(paste0('(', round(ifelse(zout[[1]]$zscore > 0, (1 - pnorm(zout[[1]]$zscore)) * 2, pnorm(zout[[1]]$zscore) * 2), 3), ')')
	, y = betas_men[1] + 0.2, x = 1, cex = 1.5)

	text(round(zout[[2]]$difference, 3), y = betas_men[2] + 0.85, x = 2, cex = 1.5)
	text(paste0('(', round(ifelse(zout[[2]]$zscore > 0, (1 - pnorm(zout[[2]]$zscore)) * 2, pnorm(zout[[2]]$zscore) * 2), 3), ')')
	, y = betas_men[2] + 0.2, x = 2, cex = 1.5)

	text(round(zout[[3]]$difference, 3), y = betas_men[3] + 0.7, x = 3, cex = 1.5)
	text(paste0('(', round(ifelse(zout[[3]]$zscore > 0, (1 - pnorm(zout[[3]]$zscore)) * 2, pnorm(zout[[3]]$zscore) * 2), 3), ')')
	, y = betas_men[3] + 0.2, x = 3, cex = 1.5)

	# Add legend
	legend('bottomleft', 
		legend = c('Women', 'Men'), 
		col = c('gray46', 'black'), 
		pch = c(17, 20), 
		lty = c(1, 1),
		bty = 'n',
		cex = 1.5)
}



genSensitivityPeriodH1 <- function(yearElection, yearPlot, file){
	pdf(file)
	outMen <- sensitivity(mydata = subset(d1, codeyear %in% centParties & ANO_ELEICAO %in% yearElection), 
		DV = 'avg_male_DF', 
		ylb = -4, 
		yub = 4, 
		dvName = 'Federal Deputy - Men', 
		clusterVar = subset(d1, codeyear %in% centParties & ANO_ELEICAO %in% yearElection)$codeyear)
	outWomen <- sensitivity(mydata = subset(d1, codeyear %in% centParties & ANO_ELEICAO %in% yearElection), 
		DV = 'avg_fem_DF', 
		ylb = -4, 
		yub = 4, 
		dvName = 'Federal Deputy - Women', 
		clusterVar = subset(d1, codeyear %in% centParties & ANO_ELEICAO %in% yearElection)$codeyear)
	zscores <- zDif(beta1 = outMen$beta, 
		beta2 = outWomen$beta, 
		se1 = outMen$se, 
		se2 = outWomen$se)$zscore
	difference <- zDif(beta1 = outMen$beta, 
		beta2 = outWomen$beta, 
		se1 = outMen$se, 
		se2 = outWomen$se)$difference
	pvalues <- ifelse(zscores > 0, (1 - pnorm(zscores)) * 2, pnorm(zscores) * 2)

	par(mar = c(5, 5, 2, 2))
	plot(x = outWomen$h[difference>0], 
		y = pvalues[difference>0], 
		pch = 20, 
		cex = 1, 
		xlim = c(0, 100),
		ylim = c(0, 1),
		xlab = 'Bandwidth', 
		ylab = 'p-value',
		axes = FALSE, 
		cex.lab = 2)
	points(x = outWomen$h[difference<=0], 
		y = pvalues[difference<=0], col = 'grey40', pch = 18) 
	segments(x0 = -5, x1 = 105, y0 = 0.05, y1 = 0.05, lty = 3)
	segments(x0 = -5, x1 = 105, y0 = 0.10, y1 = 0.10, lty = 2)
	axis(1, at = c(0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100), cex.axis = 1.5)
	axis(2, cex.axis = 1.5)
	legend('topright', 
		legend = c(expression("Men" > "Women"), expression("Men" <= "Women")), 
		col = c('black', 'grey40'),
		pch = c(20, 18), 
		bty = 'n',
		cex = 1.5)
	dev.off()
}

# This function is used to create Figure 3 in the main text, this time by party
figure3ByParty	<- function(partyName = 'PT', file){
  
  # Define Period
  d1$period <- 'weak'
  d1$period[d1$ANO_ELEICAO > 2007] <- 'strong'
  d1$period[d1$ANO_ELEICAO > 2015] <- 'moderate'
  
  # Remove other centralized party from control group
  otherCentr <- c("PC do B", "PT", "PSDB", "PPS", "PFL", "DEM", "PL", "PR", "PDT")
  otherCentr <- otherCentr[!(otherCentr %in% partyName)]
  d1 <- subset(d1, !(SIGLA_PARTIDO %in% otherCentr))
  
  
  # Define Sample
  targetParty <- d1$codeyear[(d1$SIGLA_PARTIDO %in% partyName) & d1$eleito == 1]
  
  # Weak
  out.cent.men.opt.weak <- rdrobust(y = subset(d1, codeyear %in% targetParty & period == 'weak')$avg_male_DF
                                    , x = subset(d1, codeyear %in% targetParty & period == 'weak')$run_var
                                    , p = 1
                                    , c = 0	
                                    , cluster = subset(d1, codeyear %in% targetParty & period == 'weak')$codeyear
                                    , all = TRUE)
  out.cent.women.opt.weak <- rdrobust(y = subset(d1, codeyear %in% targetParty & period == 'weak')$avg_fem_DF
                                      , x = subset(d1, codeyear %in% targetParty  & period == 'weak')$run_var
                                      , p = 1
                                      , c = 0	
                                      , cluster = subset(d1, codeyear %in% targetParty & period == 'weak')$codeyear
                                      , all = TRUE)
  
  # Strong
  out.cent.men.opt.strong <- rdrobust(y = subset(d1, codeyear %in% targetParty & period == 'strong')$avg_male_DF
                                      , x = subset(d1, codeyear %in% targetParty & period == 'strong')$run_var
                                      , p = 1
                                      , c = 0	
                                      , cluster = subset(d1, codeyear %in% targetParty & period == 'strong')$codeyear
                                      , all = TRUE)
  out.cent.women.opt.strong <- rdrobust(y = subset(d1, codeyear %in% targetParty & period == 'strong')$avg_fem_DF
                                        , x = subset(d1, codeyear %in% targetParty & period == 'strong')$run_var
                                        , p = 1
                                        , c = 0	
                                        , cluster = subset(d1, codeyear %in% targetParty & period == 'strong')$codeyear
                                        , all = TRUE)
  
  # Moderate
  out.cent.men.opt.moderate <- rdrobust(y = subset(d1, codeyear %in% targetParty & period == 'moderate')$avg_male_DF
                                        , x = subset(d1, codeyear %in% targetParty & period == 'moderate')$run_var
                                        , p = 1
                                        , c = 0	
                                        , cluster = subset(d1, codeyear %in% targetParty & period == 'moderate')$codeyear
                                        , all = TRUE)
  out.cent.women.opt.moderate <- rdrobust(y = subset(d1, codeyear %in% targetParty & period == 'moderate')$avg_fem_DF
                                          , x = subset(d1, codeyear %in% targetParty & period == 'moderate')$run_var
                                          , p = 1
                                          , c = 0	
                                          , cluster = subset(d1, codeyear %in% targetParty & period == 'moderate')$codeyear
                                          , all = TRUE)
  
  # Get point estimates, se, and ci for men
  betas_men <- c(out.cent.men.opt.weak$coef[3,], out.cent.men.opt.strong$coef[3,]
                 , out.cent.men.opt.moderate$coef[3,])
  se_men <- c(out.cent.men.opt.weak$se[3,], out.cent.men.opt.strong$se[3,]
              , out.cent.men.opt.moderate$se[3,])
  lb_men <- c(out.cent.men.opt.weak$ci[3,1], out.cent.men.opt.strong$ci[3,1]
              , out.cent.men.opt.moderate$ci[3,1])
  ub_men <- c(out.cent.men.opt.weak$ci[3,2], out.cent.men.opt.strong$ci[3,2]
              , out.cent.men.opt.moderate$ci[3,2])
  
  # Get point estimates, se, and ci for women
  betas_women <- c(out.cent.women.opt.weak$coef[3,], out.cent.women.opt.strong$coef[3,]
                   , out.cent.women.opt.moderate$coef[3,])
  se_women <- c(out.cent.women.opt.weak$se[3,], out.cent.women.opt.strong$se[3,]
                , out.cent.women.opt.moderate$se[3,])
  lb_women <- c(out.cent.women.opt.weak$ci[3,1], out.cent.women.opt.strong$ci[3,1]
                , out.cent.women.opt.moderate$ci[3,1])
  ub_women <- c(out.cent.women.opt.weak$ci[3,2], out.cent.women.opt.strong$ci[3,2]
                , out.cent.women.opt.moderate$ci[3,2])
  
  # Calculate Differences 
  zOutWeak <- zDif(beta1 = out.cent.men.opt.weak$coef[3,], beta2 = out.cent.women.opt.weak$coef[3,]
                   , se1 = out.cent.men.opt.weak$se[3,], se2 = out.cent.women.opt.weak$se[3,])
  zOutStrong <- zDif(beta1 = out.cent.men.opt.strong$coef[3,], beta2 = out.cent.women.opt.strong$coef[3,]
                     , se1 = out.cent.men.opt.strong$se[3,], se2 = out.cent.women.opt.strong$se[3,])
  zOutModerate <- zDif(beta1 = out.cent.men.opt.moderate$coef[3,], beta2 = out.cent.women.opt.moderate$coef[3,]
                       , se1 = out.cent.men.opt.moderate$se[3,], se2 = out.cent.women.opt.moderate$se[3,])
  
  # Start the plot:
  pdf(file, width = 14)
  par(mar = c(5, 5, 2, 2))
  plot(x = seq(1-0.1, 3-0.1, 1), y = betas_men, axes = FALSE, 
       pch = 20, cex = 2, xlim = c(0.75, 3.25), ylim = c(-15, 15),
       ylab = "RD Effect of Co-Partisan Mayor",
       xlab = "Electoral Rule Regime", 
       cex.lab = 2)
  sapply(1:3, function(i) lines(x = c(seq(1-0.1, 3-0.1, 1)[i], seq(1-0.1, 3-0.1, 1)[i]), 
                                y = c(lb_men[i], ub_men[i])))
  sapply(1:3, function(i) lines(x = c(seq(1+0.1, 3+0.1, 1)[i], seq(1+0.1, 3+0.1, 1)[i]), 
                                y = c(lb_women[i], ub_women[i]), col = 'gray46'))
  axis(1, at = seq(1, 3, 1), label = c('Weak\n(1998-2006)', 'Strong\n(2010-2014)', 'Moderate\n(2018)'), tick = FALSE, cex.axis = 2)
  axis(2, tick = T, cex.axis = 2, las = 2)
  segments(x0 = 0, x1 = 3.5, y0 = 0, y1 = 0, lty = 3)
  
  # Add differences
  # Weak
  segments(x0 = 0.9, x1 = 1.05, y0 = betas_men[1], y1 = betas_men[1], lty = 2)
  segments(x0 = 1.05, x1 = 1.05, y0 = betas_men[1], y1 = betas_women[1], lty = 2)
  segments(x0 = 1.05, x1 = 1.1, y0 = betas_women[1], y1 = betas_women[1], lty = 2)
  # Strong
  segments(x0 = 1.9, x1 = 2.05, y0 = betas_men[2], y1 = betas_men[2], lty = 2)
  segments(x0 = 2.05, x1 = 2.05, y0 = betas_men[2], y1 = betas_women[2], lty = 2)
  segments(x0 = 2.05, x1 = 2.1, y0 = betas_women[2], y1 = betas_women[2], lty = 2)
  # Moderate
  segments(x0 = 2.9, x1 = 3.05, y0 = betas_men[3], y1 = betas_men[3], lty = 2)
  segments(x0 = 3.05, x1 = 3.05, y0 = betas_men[3], y1 = betas_women[3], lty = 2)
  segments(x0 = 3.05, x1 = 3.1, y0 = betas_women[3], y1 = betas_women[3], lty = 2)
  
  # Add points for women's coefficients
  points(x = seq(1+0.1, 3+0.1, 1), y = betas_women, pch = 17, col = 'gray46', cex = 1.5)
  
  # Add information about difference
  text(round(zOutWeak$difference, 3), y = betas_men[1] + 2.75, x = 1, cex = 1.5)
  text(paste0('(', round(ifelse(zOutWeak$zscore > 0, (1 - pnorm(zOutWeak$zscore)) * 2, pnorm(zOutWeak$zscore) * 2), 3), ')')
       , y = betas_men[1] + 0.8, x = 1, cex = 1.5)
  
  text(round(zOutStrong$difference, 3), y = betas_men[2] + 2.45, x = 2, cex = 1.5)
  text(paste0('(', ifelse(round(ifelse(zOutStrong$zscore > 0, (1 - pnorm(zOutStrong$zscore)) * 2, pnorm(zOutStrong$zscore) * 2), 3) == 0, "0.000", round(ifelse(zOutStrong$zscore > 0, (1 - pnorm(zOutStrong$zscore)) * 2, pnorm(zOutStrong$zscore) * 2), 3)), ')')
       , y = betas_men[2] + 1.05, x = 2, cex = 1.5)
  
  text(round(zOutModerate$difference, 3), y = betas_men[3] + 2.2, x = 3, cex = 1.5)
  text(paste0('(', round(ifelse(zOutModerate$zscore > 0, (1 - pnorm(zOutModerate$zscore)) * 2, pnorm(zOutModerate$zscore) * 2), 3), ')')
       , y = betas_men[3] + 0.8, x = 3, cex = 1.5)
  
  # Add legend
  legend(x = 0.7, y = -7.5, 
         legend = c('Women', 'Men'), 
         col = c('gray46', 'black'), 
         pch = c(17, 20), 
         lty = c(1, 1),
         bty = 'n',
         cex = 1.5)
  
  dev.off()
}


